oy^^c!t,  OF,  Canada 


ASfiOMAUTlCAL  REPORT 


•  -■*>  i.  ■  ^■>g  "  ■  ..J>'  '’y^-'  Trk  '.  ':  V'-i^vC'' »'■ »‘  '  •  *'.  '  ►.V:  .;'  ^ • 


CLOSED  FORM,  FINITE  ELEMENT  SOLUTIONS  FOR  PLATE  VIBRATIONS 


by 

GARRY  M.  LINOBERG,  MERVYN  0.  OLSON  AND  HELEN  A.  TULLOCH 


A,  H,  Hall,  Head 

Structures  and  Materials  Section 


F.  K.  Tburston 
Director 


SUMMARY 


The  convergence  rates  of  eigenvalue  solutions  using 
two  finite  plate-bending  elements  are  studied.  The  elements 
considered  are  the  well-known  12-degree-of-freedom,  non- 
conforming  rectangular  element  and  the  16-degree-of-freedom, 
conforming  rectangular  element.  Three  problems  are  analyzed: 
a  square  plate  simply  supported  on  two  opposite  sides,  with 
the  other  two  sides  clamped,  simply  supported,  or  free.  Closed 
form,  finite  element  solutions  for  these  problems  are  obtained 
by  using  shifting  E-operators. 

With  few  exceptions,  eigenvalue  solutions  found 
with  the  non-conforming  element  converge  from  below  the 
exact  answers  at  an  asymptotic  rate  of  n”",  where  n  is  the  number 
of  elements  on  a  side.  However,  since  the  array  size  needed  for 
such  convergence  is  very  large,  little  can  be  said  about  the 
convergence  rates  for  practical  arrays.  The  conforming  element 
solutions  converge  from  above  at  a  rate  of  n  for  values  of  n 
larger  than  6.  A  comparison  of  the  errors  involved  in  using 
these  two  elements  shows  that  the  conforming  element  is  far 
superior  to  the  non-conforming  element. 
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CLOSED  FORM,  FINITE  ELEMENT  SOLUTIONS  FOR  PLATE  VIBRATIONS 
1.0  INTRODUCTION 

In  recent  years,  a  new,  powerful  approximate  technique  for  solving  complicated  boundary 
value  problems  has  evolved.  This  is  the  so-called  finite  element  method,  and  it  has  been  extensively 
developed  for  the  analysis  of  structures.  In  this  finite  element  technique,  the  continuous  system  is 
replaced  by  a  substitute  system  consisting  of  a  number  of  finite  elements  linked  together.  Once  the 
properties,  stiffness,  mass,  etc.  of  the  individual  elements  have  been  defined,  the  whole  substitute 
system  can  be  described  by  large  matrix  equations,  readily  solvable  using  modern  computers. 

The  success  of  the  finite  element  technique  depends  upon  the  behaviour  of  each  element 
and  on  how  few  elements  it  takes  to  adequately  model  a  real  structure.  It  is  important  that  the  finite 
element  solutions  converge  to  the  exact  answers  as  the  number  of  elements  is  increased,  and  that  thus 
convergence  be  rapid.  It  is  also  desirable  to  have  some  estimate  of  the  accuracy  of  any  solutions 
found  using  a  specific  array  of  elements. 

For  plate-bending  problems,  two  of  the  most  important  rectangular  elements  developed 
are  the  12-degree-of-freedom  non-conforming  element  and  the  16-degree-of-freedom  conforming 
element.  The  first  element,  independently  derived  by  several  investigators,  Lindberg”,  Melosh^\ 
Zienkiewicz  and  Cheung'”,  Dawe'",  and  Clough  and  Tocher’”  amongst  others,  was  one  of  the  first 
successful  plate-bending  finite  elements  developed.  This  element  is  called  non-conforming  since  two 
elements  linked  together  will  have  continuous  displacements  along  their  common  edge,  but  will  not 
have  continuous  normal  slopes.  More  recently,  a  16-degree-of-freedom  element  has  been  inde¬ 
pendently  derived  by  Bogner,  Fox  and  Schmit'”,  Butlin  and  Leckie’\  and  later  by  Mason”’.  This 
element  has  the  added  property  that  two  elements  linked  together  have  continuous  normal  slopes 
along  their  common  edge  as  well  as  continuous  displacements. 

Recent  literature  has  clarified  to  some  extent  the  criteria  that  assure  convergence  of  the 
finite  element  approximation  to  the  true  solution  as  the  array  is  refined.  A  sufficient  condition  is 
that  the  element  be  capable  of  representing  a  constant  strain  condition;  in  the  case  of  plate  bending 
this  should  be  a  constant  curvature  condition,  pure  bending,  or  pure  twist.  This  includes  a  condition 
of  zero  strain  that  corresponds  to  rigid  body  displacements.  A  further  condition  is  that  the  element 
be  conforming.  This  condition  assures  monotonic  convergence  of  the  total  potential  energy  (Ref. 
9,10,11).  Both  elements  considered  in  this  study  satisfy  the  first  criterion  and,  hence,  yield  results 
that  must  converge  to  correct  solutions.  The  16-degree-of-freedom  element  also  satisfies  the  second 
criterion  and,  hence,  must  provide  monotonic  convergence  of  potential  energy  for  static  problems. 
Indeed,  it  is  possible,  following  the  arguments  presented  by  Cowper  et  al."’,  to  show  that  the  rate  of 
convergence  should  be  proportional  to  n""*,  where  n  is  the  number  of  elements  on  the  side  of  an 
array.  However,  at  this  time,  an  extension  of  this  general  convergence  proof  to  cover  dynamic 
problems  is  still  lacking. 

The  convergence  of  the  12-degree-of-freedom  model  has  be<  n  studied  by  Walz,  Fulton  and 
Cyrus’'^’  who  show  that  for  simply-supported  square  plates,  both  the  static  and  dynamic  conver¬ 
gence  rates  should  be  n"l  For  large  values  of  n,  they  also  provide  an  estimate  of  the  errors  involved 
in  using  these  elements.  This  convergence  study  is  somewhat  limited,  since  only  one  boundary 
condition  is  considered,  and  the  results  hold  true  only  for  large  arrays. 
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A  more  complete  study  of  the  dynamic  convergence  rates  of  these  two  elements  is  given  in 
this  report.  Three  problems  are  studied;  a  square  plate  simply  supported  on  two  opposite  sides, 
with  the  other  two  sides  simply  supported,  clamped,  and  free.  The  problems  chosen  all  have  exact 
solutions,  so  that  the  finite  element  errors  may  be  systematically  studied.  A  closed  form  type  of 
solution  using  finite  shifting  E-operators  is  developed,  so  that  solutions  for  large  arrays  can  be 
found  with  little  computational  effort.  This  type  of  solution  also  reveals  results  of  general  interest 
that  might  easily  be  masked  by  a  mass  of  arithmetical  detail.  Twenty  eigenvalues  are  found  for 
each  of  the  problems,  using  each  of  the  elements.  The  number  of  elements  on  a  side  is  varied  from 
2  to  20,  and  the  effects  of  reducing  boundary  conditions  on  some  edges  are  investigated. 

2.0  EXACT  SOLUTIONS 

The  method  for  finding  exact  solutions  to  these  three  problems  will  now  be  given.  The 
well-known  differential  equation  governing  the  free  vibrations  of  plates  is  (Timoshenko'  ”,  p.  334) 

vw  -  - 

f^t' 

where  W(x,y,t)  is  the  deilection  of  the  plate,  /j  is  the  mass  per  unit  area  of  the  plate,  D  is  the  plate¬ 
bending  rigidity,  and  t  is  the  time.  Assuming  simple  harmonic  motion  of  the  plate 

W  (x,  y,  t)  =  w  (x,  y)  sin  wt  (2) 

and  substituting  this  into  equation  (1)  leads  to 

^  4.  0  4.  ^  w 

<lx*  "  dx^dy'  0y^  D 

Consider  the  plate  shown  in  Figure  1,  with  the  two  opposite  sides  y  =  0  and  y  =  L  simply 
supported.  These  boundary  conditions  are  satisfied  by  a  sine  wave  in  the  y-direction,  so  that 

w  (X,  y)  =  X  (x)  sin  (p7ry/L)  (4) 

where  p  is  the  number  of  half  waves  in  the  y  direction.  Thus,  the  differential  equation  to  be  solved 
becomes 


X""  -  2  X"  +  X  =  0  (5) 

where  X'"'  =  0'X/^x^  etc. 

Thus,  the  partial  differential  equation  has  been  reduced  to  a  solvable  ordinary  differential 
equation.  Assuming  that 

X  (x)  =  A  exp  («7rx/L)  (6) 

and  substituting  this  into  equation  (5)  yields  the  characteristic  equation 


a'  -  2p’a*  +p*  (1  -  «’)  =  0 


(7) 


wliere  w'  =  >ita'-'LVpV'D  is  a 

non-dimensional  frequency  parameter. 

Solving  for  the  roots  of  this 

(luartic  ecjuntion  yields 

If  u  >  1,  then 

a*  =  p*  (1  i  w) 

(8) 

rt  =  ±  k  ,  where  k*  =  p“  (1  -f-  5) 

(9) 

and 

a  =  ±  im  ,  where  m*  =  p’  (w  —  1) 

(10) 

and  i  is  the  imaginary  unit. 

Using  these  four  roots,  the  solution  for  X  becomes 

X(x)  =  A  cosh(kjrx/L)  +  B  sinh(k7rx/L)  +  C  cosCm^x/L)  +  E  8in(mirx/L)  (11) 

If  w  <  1,  then 

«  =  ±  k  ,  where  k*  =  p’  (1  +  w)  (12) 

and  a  =  ±  where  =  p*  (1  -  C)  (13) 

so  that 

X(x)  =  F  cosh(kirx/L)  +  G  sinh(k7rx/L)  +  H  cosh  (/ttx/L)  +  J  sinh(/7rx/L)  (14) 

The  constants  A,  B,  C,  and  E  or  F,  G,  H,  and  J  can  be  determined  by  consider¬ 
ation  of  the  boundary  conditions  along  the  two  sides,  x  =  ±L/2.  Two  boundary  conditions  on 
each  side  are  used  to  give  four  equations  in  the  four  unknown  constants,  and  elimination  of  these 
gives  a  frequency  equation  for  the  particular  case. 

It  is  easier  to  consider  symmetric  and  anti-symmetric  solutions  separately.  For  solutions 
that  are  symmetric  in  x  about  the  y-axis,  equation  (4)  becomes 

w(x,y)  =  [A  cosh(k7rx/L)  C  cos(m7rx/L)]  sin(p7ry/L)  (15) 

while  for  the  anti-symmetric  case  it  becomes 

w(x,y)  =  IB  sinh(k7rx/L)  +  E  sin(m7rx/L)]  sin(p7ry/L)  (16) 

Boundary  conditions  for  the  three  types  of  edges  considered  are  well  known.  For  the  clamped 
edge,  the  boundary  conditions  are 

w  (x,  y)  j  =0 

I  x=L/2 


dW  =  0 

ax  X  =  L/2 


and 


(17) 
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for  the  simply-supported  edge,  they  are 


w  (X,  y)  =0 

X  =  L/2 


and  M. 

Vax-  ay/ 

and  for  the  free  edge,  they  are 

=  + 

\dx'  dy/ 


=  0 


X  =L/2 


M* 


=  0 


and 


X  =L/2 

Q  .  (v  -  .  (i!!?  +  ,2  -  ,) 

\  I  \dx  axdy- 1 


=  0 


X  =L/2 


(18) 


(19) 


The  frequency  equations  for  three  types  of  boundary  conditions  are  obtained  by  substituting 
equation  (15)  into  the  two  boundary  conditions,  and  then  eliminating  the  constants  from  the  ensuing 
equations.  The  transcendental  equations  are 


and 


k  tan(m7r/2)  -f  m  tanh(k7r/2)  =  0 
cosh(k7r/2)  cos(mir/2)  =  0 

m  tan(m7r/2)  +  k  tanh(k7r/2)  =  0 

\  +  I  \  k*  -  ^p-  / 


(20) 

(21) 

(22) 


for  the  clamped,  simply-supported,  and  free  boundaries,  respectively. 

Similar  transcendental  equations  can  be  found  for  the  anti-symmetric  cases  and  for  the 
cases  where  w  <  1. 


The  eigenvalues  for  the  problems  correspond  to  the  roots  or  zero  values  of  the  trans¬ 
cendental  equations,  and  may  be  found  by  an  iterative  process.  A  value  of  p,  the  number  of  half 
waves  in  the  y-direction,  is  selected  and  values  of  the  non-dimensional  frequency  parameter  w 
are  used  to  evaluate  the  transcendental  equation  being  solved.  When  a  zero  crossing  is  found,  an 
iterative  procedure  is  used  to  obtain  a  precise  value  of  frequency  that  makes  the  equation  as  close 
to  zero  as  required.  The  exact  non-dimensional  frequencies  found  for  these  three  problems  are 
given  later. 


3.0  CLOSED  FORM  FINITE  ELEMENT  SOLUTIONS 


The  straightforward  method  of  obtaining  eigenvalues  for  various  finite  element  gridworks 
is  to  set  up  the  large  system  of  simultaneous  equations  for  a  particular  gridwork,  apply  boundary 
conditions,  and  solve  the  resulting  eigenvalue  problem  on  a  digital  computer.  There  are  two  diffi¬ 
culties  in  using  such  a  procedure  for  the  present  study.  Firstly,  the  size  of  the  eigenvalue  problem 
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risos  rapidly  with  increase  of  the  number  of  elements  used,  and  hence  the  present  study  would  in¬ 
volve  a  prohibitive  amount  of  computation  time.  Secondly,  results  of  general  interest  can  be  hidden 
in  a  mass  of  arithmetic  detail. 

A  method  of  closed  form  solution  utilizing  shifting  E-operators  has  therefore  been  adopted 
for  this  study.  This  method  was  developed  by  Lindberg"  in  a  comparison  study  of  several  different 
finite  elements.  It  was  also  used  in  a  simplified  form  by  Leckie'^'.  The  procedures  are  similar  to  those 
used  in  finding  the  exact  solutions  for  the  three  problems.  The  equilibrium  equations  for  an  internal 
point  of  an  as.semblage  of  elements  can  be  written  using  the  stiffness  and  mass  matrices  of  the  ad¬ 
jacent  elements.  These  equilibrium  equations  are  expressed  in  terms  of  the  generalized  displacements 
of  the  surrounding  element  corners,  and  shifting  E-operators  are  used  to  express  these  equations  in 
terms  of  the  generalized  displacements  of  the  single  internal  point.  Expressions  for  these  generalized 
displacements,  sinusoidal  in  the  direction  normal  to  the  simple  supports  and  exponential  in  the  other 
direction,  are  assumed  and  substituted  into  the  equilibrium  equations.  A  characteristic  equation 
for  the  assumed  displacement  functions  is  then  obtained,  and  general  expressions  for  the  displace¬ 
ment  functions  are  found.  These  general  functions  can  then  be  substituted  into  the  generalized 
force  equations  for  an  edge  point,  thus  obtaining  both  force  and  displacement  equations  for  an  edge. 
By  selecting  the  appropriate  boundary  conditions,  an  approximate  transcendental  frequency  equation 
for  each  of  the  three  cases  is  found.  The  roots  of  these  equations  are  found  using  iteration  procedures. 

These  equations  are  a  function  of  the  number  of  elements  used  in  the  problem,  and  it  is  easy 
to  vary  this  number  and  study  the  convergence  of  the  approximate  solutions.  It  is  also  possible  to  use 
a  large  number  of  elements  without  increasing  the  amount  of  computing  time  required  for  a  solution. 

3.1  Twelve-degree-of-freadom,  Non*conforming  Element 

This  element  has  three  degrees  of  freedom  per  corner,  ^y,  and  w/a,  a  total  of  12  degrees 
of  freedom.  It  is  commonly  called  a  non-conforming  element,  since  two  elements  linked  together 
have  continuous  displacements  along  their  common  edge,  but  do  not  have  continuous  normal  slopes 
there.  The  element  is  well  described  in  the  literature  (Ref.  1-5),  so  only  the  numerical  values  of  the 
stiffness  and  mass  matrices  used  in  this  study  need  be  given.  These  stiffness  and  mass  matrices  for 
a  square  element  of  side  a,  and  Poisson’s  ratio  of  are  given  in  Tables  1  and  2  respectively,  where 
the  displacement  vector  is 

{4'xi,  fvi.  wj/a,  . W4/a) 

The  problem  to  be  solved  is  shown  in  Figure  2  for  the  particular  element  assemblage  of 
n  =  2q  =  6.  Note  that  the  co-ordinates  r  and  s  are  not  continuous,  but  rather  are  incremental  and 
refer  to  assemblage  points.  Consider  the  internal  point  5  of  the  assemblage  shown  in  Figure  3.  Since 
the  displacements  at  the  internal  point  are  assumed  to  be  continuous,  it  is  possible  to  write  a  matrix 
equation  for  the  forces  at  this  point  as: 

17  0  39  44  0  0  17  0  -39  56  0  192  272  0  0 

0  17  39  0  56  192  0  17  39  0  44  0  0  272  0 

_39  _  39  -  66  0  -192  -  408  39  -  39  -  66  -192  0  -  408  0  0  1896 


56 

0 

-192 

17 

0 

39 

44 

0 

0 

17 

0 

-39 

0 

44 

0 

0 

17 

-  39 

0 

56 

-192 

0 

17 

-39 

(W, 

192 

0 

-408 

-39 

39 
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192 

-408 

39 

39 

-66 
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nu*&* 

25200 


-30  -  28-116 
-28  -  30-116 
116  116  394 


80  0  0  -30  28  116  -120  0  -548 

0  -120  -548  28  -30  -116  0  80  0 

0  548  2452  -116  116  394  548  0  2452 


320 

0 

0 


0 

320 

0 


0 

0 

13816 


where 


-120 

0 

548 

-30 

28 

-116 

80 

0 

0 

-30 

-28 

116 

0 

80 

0 

28 

-30 

116 

0 

-120 

548 

-28 

-30 

116 

-548 

0 

2452 

116 

-116 
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0 

-548 

2452 
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-116 
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jW,)'*'  =  Wi/a.  ^r2,  ^1^,2,  W2/a, . ^p,o,  Wg/a} 


(W,)  (23) 


is  a  27-component  vector  of  the  displacements  at  assemblage  points  surrounding  point  5. 


Since  only  dynamic  solutions  with  no  external  loads  are  required,  then  for  equilibrium  at 
any  point  these  forces  must  equal  zero. 

A  shifting  E-operator  (see  App.  A)  is  defined  as 

Ejf  Er  w,.,  =  w,4,k.,+„, 

where  w,^  is  the  displacement  w  at  the  point  (r,s)  in  Figure  3.  Then,  using  the  r  and  s  co-ordinates 
for  points  1  to  9  shown  in  Figure  3,  these  equilibrium  equations  may  be  written  as 

{  (Er‘  +  E‘)  (E.-‘  -f  Ei)  (17  -f-  30  7*)  -f  E?  (E,"*  +  E\)  (44  -  80  7*) 

-f  E2  (Er‘  +  Ei)  (56  -f  120  7*)  +  E?  E°  (272  -  320  7*)  } 

-f  I  (Er‘  -  Ei)  (Er‘  -  Ei)  (28  7*)  }  ^.5 

-I-  {  (Er‘  +  Ei)  (E7‘  +  Ei)  (39  +  116  7*)  +  E”  (Er‘  -  Ei)  (192  -f  548  7*)  }wr./a  =  0  (24) 

{  (Er‘  -  Ei)  (E.-'  -  Ei)  (28  7*)  1  ^,5 

-f  I  (Er‘  4-  Ei)  (E.-‘  +  Ei)  (17  4-  30  7*)  +  E?  (E,"'  4-  Ei)  (56  +  120  7*) 

4-  E2(Er‘  4-  Ei)  (44  -  80  7*)  +  E?  E°  (272  -  320  7*)  !  >^.5 

+  1  (Er‘  4-  Ei)  (E.-‘  -  Ei)  (39  f  116  7*)  4-  E?  (E,"’  -  Ei)  (192  4-  548  7*)  }  w.,/,  =  0  (25) 

(  (E7‘  -  Ei)  (E.-'  4-  Ei)  (-  39  -  116  7*)  -  E®(Er‘  -  Ei)  (192  4-  548  7*)  1  lArr. 

4-  {  (Er’  4-  Ei)  (E.-‘  -  Ei)  ( -  39  -  116  7*)  -  E?(Er‘  -  Ei)  (192  -h  548  7*)  1  ^.5 

4-  {  (Er‘  4-  Ei)  (E.-‘  4-  Ei)  ( -  66  -  394  7*)  -  E2(Er‘  4-  Ei)  (408  4-  2452  7*) 

-  E?(E,-'  4-  Ei)  (408  2452  7*)  +  E"  Ej  (1896  -  13816  7*)  1  Wg/,  =  0 


In  these  expressions,  7* 


=  liu  a  /  560  D  is  a  non-dimensional  frequency  parameter. 


(26) 
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Since  the  problems  considered  are  all  simply  supported  on  two  opposite  sides,  it  is  possible, 
as  in  the  exact  solution,  to  assume  that  the  deflected  shape  in  the  s-direction  is  a  sine  wave.  If  it  is 
assumed  that  the  deflection  in  the  r-direction  varies  exponentially,  then 

(w/a)r.,  =  A  e^'  sin(p7r8/n) 

(^r)r.i  =  B  e"'  sin(p7r8/n)  (27) 

=  C  e^''  co8(pir8/n) 

where  p  is  the  number  of  half  sine  waves  in  the  s-direction. 

One  of  the  valuable  properties  of  E-operators  is  that  they  obey  the  rule 
Fi(Er)  Fa  (E.)  a*''  b™  =  a**'  b"“  Fi  (a‘‘)  Fa  (b™) 

where  Fj,  Fa  are  specified  functions.  Hence,  it  is  easy  to  substitute  exponential  functions  into 
E-operator  equations.  A  table  of  E-operator  transformations  is  given  in  Appendix  A.  Substituting 
the  deflection  equations  (27)  into  the  equilibrium  equations  (24-26)  and  using  this  transformation 
table  yields  the  following  three  equations 

sinh  <r  1  -39  cos(pir/n)  -  96  -  7  (58  cos(p7r/n)  -4-  137)  ]  A 

-f  (cosh  c  [17  cos(p7r/n)  -f-  28  -f  7  (15  cos(pir/n)  +  30)  ) 

-f  22  co8(pir/n)  -f  68  —  7  (20  co8(prr/n)  40))  B 

—  14  7  sinh  «r  sin(p7r/n)  C  =  0  (28) 

8in(pir/n)  (cosh  <r  ( —39  —  687)  —  96  —  137  7  )  A  -f  14  7  sinh  a  sin(p7r/n)  B 
+  (cosh  <T  [17  cos(p7r/n)  -f  22  -4-  7  (15  cos(p7r/n)  -  20)]  +  28  cos(p7r/n)  -|-  68 
-h  7  (30  cos(pir/n)  —  40)  }  C  =  0  (29) 

(cosh  a  1  —66  co8(p7r/n)  —  204  —  7  (197  cos(p7r/n)  +  613)  ] 

—  204  co8(pir/n)  -f  474  —  7  (613  cos(p7r/n)  -4-  1727)  }  A 

-4-  (sinh  a  [39  C08(pir/n)  -f-  96  +  7  (58  C08(pir/n)  -4-  137)]  )  B 

-f-  8in(pir/n)  (cosh  <r(  —  39  —  887)  — 96  —  137  7  }  C  =  0  (30) 


where  7  =■  y*/2. 


—  8  — 


For  a  nontrivial  solution  for  A,  B,  and  C,  the  determinant  of  the  coefficients  must  vanish, 
and  this  determinant  simplifies  to  an  equation  cubic  in  cosh  a.  This  is  the  characteristic  equation 
for  these  problems,  and  its  three  roots  give  three  values  for  a  that  can  be  used  to  determine  w.  For 
these  problems,  the  roots  are  found  to  be 

cosh  ^  1 
—  1  <  cosh  tr  ^  1 

\ 

and  cosh  a  <  —  1 

If  cosh  ffi>  1,  then  ai  =  ±  «  say,  and  so 


I  (w/a),.„l,  =  (E  cosh  ar  +  F  sinh  «r)  sin(p7rs/n)  (31) 

If  —  1  ^  cosh  a-i  ^  1,  then  n-,  =  +  i^  say,  and  so 

Kw/aV.*],  =  (G  cos  Hr  +  H  sin  /3r)  sin(p7rs/n)  (32) 

Lastly,  consider  cosh  tr.-i  <  —  1.  Say  that  cosh  03  =  -6,  so  that  cru  =  cosh“*(S)  +  iff  =  {  +  iff.  Then 
(  (w/a)r.,  I,  =  (I  (  -1)'  cosh  {r  4-  J  (-!)'•  sinh  {r)  sin(pffs/n)  (33) 

Adding  the  three  components  gives  a  final  expression  for  the  deflection  as 
(w/a)r,*  =  {  E  cosh  or  +  F  sinh  or  +  G  cos  )3r  +  H  sin 

+  I  ( —  !)'■  cosh  ^r  +  J  f  —  !)'■  sinh  {r  }  sin(pffs/n)  (34) 


It  is  interesting  to  note  that  the  first  four  terms  are  similar  to  those  found  in  the  exact  solution, 
while  the  last  two  terms  arise  as  a  result  of  the  finite  element  approximation. 

As  in  the  exact  solution,  it  is  easiest  to  consider  symmetric  and  anti-symmetric  cases 
separately.  If  only  symmetric  solutions  are  considered,  the  deflection  becomes 

(w/a)r..  =  (  E  cosh  or  -f  G  cos  /3r  -f  I  ( -  I)*"  cosh  {r  }  sin(pffs/n)  (35) 

By  using  the  first  two  equilibrium  equations,  it  is  possible  to  solve  for  (^,),,,  and  (i^,),,,  in 
terms  of  (w/a),,.  If  the  above  expression  for  the  deflection  is  substituted  into  these  expressions,  and 
E-operator  transformations  are  applied,  then  these  expressions  become 

(^r)r.»  =  (  E  Xi  sinh  or  -h  G  X2  sin  /3r  -f  I  Xa  ( -  I)*"  sinh  {r  }  sin(pffs/n)  (36) 

(^.)r.»  =  {  E  X<  cosh  or  -f  G  Xr,  cos  0r  -i-  I  Xo  ( -  1)'  cosh  {r  }  cos  (pffs/n)  (37) 

where  the  X  terms  are  complicated  functions  of  a,  i3,  {  and  ^ .  These  terms  are  given  in  Appendix  B. 
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Now  the  boundary  conditions  must  be  considered.  The  finite  element  assemblage  for  an 
edge  point  is  shown  in  Figure  4.  Using  the  stiffness  and  mass  matrices  of  these  ejements,  it  is  possible 
to  write  the  equilibrium  equations  for  this  edge  point  as 


)  p. 
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56 
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-120 
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28  -116 
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0 

160  0 

28 

-30  116 
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(38) 

548 

0  2452  - 

922 

0  6908 

116 

-116  394 

-199  -274  1226 

where 

{W2)  =  {\Ari,  Wi/a,  \P,2,  ^.2,  W2/a, . ,  V'rfi,  We/a} 


Using  E-operators,  these  equations  for  and  (V»a)„.,  can  be  written  in 

terms  of  (^J,,,„,and  (w/a),,,H.  The  expressions  for  these  edge  displacements  are,  since  r  =»  q  -  n/2, 

=  (E  Xi  sinh  aq  -t-  G  X2  sin  /3q  4-  I  X3  ( —  1)''  sinh  {q}  sin(p7rs/n) 

==  {EX4  cosh  aq  +  G  X.-,  cos  <3q  +  I  Xr  ( —  1)''  cosh  {q}  cos  (p7rs/n)  (39) 

(w/a)„,«  =  {E  cosh  oq  +  G  cos  0q  +  I  ( -1)"  cosh  fq)  sin(p7rs/n) 

Substituting  these  into  the  expressions  for  edge  forces  and  using  E-operator  transformations  gives 

(Mr),,.,.  =  (<i  E  +  €2  G  -I-  1}  sin(p7rs/n) 

tMn),,.,  =  ~{«4  E  +  G  +  ««  I)  cosCpTTs/n)  (40) 

(V*a)„.,  =  ^(«7  E  -f-  €8  G  +  ej)  I)  sin(p7rs/n) 

where  the  «’s  are  complicated  functions  of  «,  ^  and  aq,  t)ci,  {q.  and  pir/n.  These  functions  are  also 

given  in  Appendix  B. 

Now  both  the  edge  forces  and  the  edge  displacements  have  been  defined  (eq.  (39)  and  (40)  ). 
By  selecting  the  appropriate  boundary  conditions  from  these,  transcendental  equations  can  be 
found  that  yield  the  required  eigenvalues  of  the  three  problems.  These  boundary  conditions  will 
now  be  considered. 
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Since  there  are  tliree  unknowns  in  the  edge  expressions,  three  boundary  conditions  must 
he  satisfied.  This  is  different  than  the  exact  solutions  where  only  two  boundary  conditions  could  be 
satisfied,  and  arises  because  this  type  of  solution  is  approximate.  In  the  continuous  case,  if  a  function 
such  as  the  displacement  w  is  set  equal  to  zero  along  an  edge,  then  all  the  derivatives  of  w  taken 
tangent  to  that  edge  are  also  automatically  set  equal  to  zero.  This  does  not  automatically  follow  in 
approximate  solutions  and,  hence,  it  is  reasonable  to  have  to  satisfy  additional  boundaiy  conditions. 

The  boundary  conditions  for  a  clamped  edge  are 

(w/a).,,«  =  0  ;  =  (0  ;  =0  (41) 

Note  that  in  the  continuous  case  the  third  boundary  condition  (^,),,,h  =  0  would  automatically  be 
satisfied  by  prescribing  the  fn-st  boundary  condition,  (w/a\,,„.=  0,  while  in  the  approximate  case 
this  is  not  true. 

For  a  simply-supported  edge,  the  boundary  conditions  are 

(w/a>„,g  =  0  ;  (Mr),,.*  =  0  ;  (i/'*).,,*  =  0  (42) 

Again,  the  first  and  third  boundary  conditions  are  intimately  related.  Finally,  for  a  free 
edge,  the  boundary  conditions  are 

(Mr),,.*  =  0  ;  (M*)„.,  =  0  ;  (V*a).,.,  =  0  (43) 

The  difference  between  the  boundary  conditions  for  the  continuous  case  and  those  for  the  approxi¬ 
mate  case  is  greatest  for  this  last  problem.  In  the  continuous  case,  a  compound  shear  condition  was 
set  equal  to  zero.  This  consisted  of  the  shear  force  plus  the  rate  of  change  of  twisting  moment.  Here, 
both  the  shear  force  and  the  tangential  bending  moment  are  set  equal  to  zero  (the  second  and  third 
conditions)  but  this  tangential  bending  moment  does  not  correspond  to  the  twisting  moment  of  the 
continuous  theory. 

For  each  problem,  the  appropriate  edge  expressions  are  set  equal  to  zero,  and  the  con¬ 
stants  E,  G,  and  I  are  eliminated  to  yield  a  transcendental  frequency  equation.  As  for  the  exact 
solutions,  iteration  procedures  are  used  to  find  the  roots  of  these  expressions,  and  hence  the  eigen¬ 
values.  Note  that  each  transcendental  equation  is  a  function  of  q,  the  number  of  finite  elements  on 
a  half  side  of  the  problem.  Thus,  eigenvalues  can  be  found  for  any  desired  value  of  q,  and  hence  n, 
so  it  is  possible  to  study  the  convergence  of  solutions  as  n  is  varied. 


3.1.1  Two-root  Approximoto  Solutions 

In  the  exact  solutions,  only  two  boundary  conditions  on  an  edge  may  be  .satisfied.  It  has 
been  shown  that  three  conditions  must  be  satisfied  for  the  approximate  solutions.  An  investigation 
was  carried  out  to  find  the  effect  of  leaving  out  the  third  boundary  condition. 

Consider  the  deflection  expression,  equation  (39) 

(w/a)r.*  =  IE  cosh  ar  -f-  G  cos  j3r  +  I  (  —  1)'  cosh  {r]  sin(p7rs/n) 
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1  he  fij'st  (wo  terms  are  similar  lo  those  used  in  the  exact  solution,  while  the  third  term  appears  to 
he  different.  Assuming  that  it  can  be  disregarded,  the  deflection  becomes 

(w/a)r.«  =  IE  cosh  fltr  +  G  cos  ,flr]  8in(p7rs/n)  (44) 

and  slopes  and  edge  forces  are  the  same  as  before,  with  the  third  term  dropped.  The  boundary  con¬ 
ditions  used  in  the  two-root  solutions  are; 

for  the  clamped  edge 


(w/a),,.,  =  0  ;  (^r).,.,  =  0 

(45) 

for  the  simply-supported  edge 

(w/a),,.,  =  0  ;  (Mr),,.,  =  0 

(46) 

and  for  free  edge 

0;i(V.a)„.+  .0 

(47) 

Note  that,  in  the  free  edge,  the  boundary  condition  used  only  approximates  the  classical  boundary 
condition.  As  before,  the  application  of  the  boundary  conditions  yields  transcendental  equations. 
The  solutions  of  these  two-root  cases  will  be  discussed  later. 

3.2  Sixtean-degree-of-fresdom  Conforming  Elemont 

This  element,  independently  derived  by  Bogner,  Fox  and  Schmit*^’,  Butli\<  and  Leckie'*, 
and  Mason'",  is  called  conforming  since  elements  linked  together  have  continuc'.'s  normal  slopes 
along  their  common  edge  as  well  as  continuous  displacements.  There  are  four  degrees  ov  freedom  at 
each  corner  point,  'Py,  a^^y,  and  w/a.  The  stiffness  and  mass  matrices  for  a  square  element  of 
side  a,  and  Poisson’s  ratio  of  li  are  given  in  Tables  3  and  4,  where  the  displacement  vector  is 

{wi/a,  iyi,  a  ^,yi,  w^/a, . ,  a  \p^y^} 


Closed  form  solutions  using  this  element  are  foucid,  <'ollowing  the  same  procedure  as 
before.  Four  equilibrium  equations  for  an  internal  point  are  written,  and  by  assuming  displacement 
functions  similar  to  equations  (27),  i.e. 

(w/a)r,*  =  A  e"’’  sin(p7rs/n),  etc. 

and  applying  E-operator  transformations,  a  characteristic  equation  for  the  problems  is  found.  This 
equation  is  quartic  rather  than  cubic  as  before,  and  its  four  roots  give  values  for  a  that  are  used  to 
determine  (w/a).  For  these  problems,  it  is  found  that  all  four  roots  are  greater  than  —1.  Typically 

cosh  ffi  ^  1  ;  (Ti  —  i  a 

cosh  ff?  ^  1  ;  (72  =  ±5 

cosh  (73  ^  1  ;  <T.i  =  +  <^» 

- 1  '.rS  cosh  <74  ^  1  ;  (74  =  ±  i/3 
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which  yields  the  deflection  expression 

(w/a)r,,  =  I A  cosh  ar  +  B  sinh  «r  +  C  cosh  5r  +  D  sinh  5r  +  E  cosh  tt>r 
4-  F  sinh  <f>v  +  G  cos  /3r  4-  H  sin  |3i'|  sin(p7rs/n) 


(48) 


In  some  cases,  two  of  the  roots  are  between  -1  and  1,  and  the  deflection  expression  then  contains 
four  hyperbolic  expressions  and  four  trigonometric  expressions. 


As  before,  expressions  for  the  slopes  and  twist  are  found  in  terms  of  the  deflection  using 
the  equilibrium  equations.  The  expression  for  (w/a)  is  substituted  into  these  equations  and,  asing 
E-operator  transformations,  they  become 

(^r)r.«  =  lA  Xio  sinh  ar  4-  C  Xn  sinh  5r  4-  E  X12  sinh  ii>r  4-  G  Xm  sin  0r]  sin(p7rs/n) 

(^»)r.ii  =  [A  Xi4  cosh  ar  4-  C  Xi.*-,  cosh  Sr  4-  E  Xm  cosh  ^r  4-  G  X17  cos  lir\  cos(p7rs/nl  (49) 

(a^r«)r..  =  [A  Xis  sinh  ar  4-  C  Xis»  sinh  Sr  4-  E  X20  sinh  (j>r  4-  G  X21  sin  iir\  cos(p7rs/n) 


where  only  symmetric  terms  are  considered.  No  details  of  tlie  equations  or  definition  of  the  X’s  are 
given,  since  the  expressions  involved  are  exct'edingly  cumberso  me.* 

Boundary  conditions  are  then  considered  as  before,  and  expressions  for  the  edge  forces 
are  found  as 


(Mr)(|,i.  =  D  [tio  A  4"  €ji  C  4-  €12  E  4-  *13  G)  sin(p7rs/n) 

(M,)q,„  =  D  lei4  A  4-  €|5  C  4-  ei6  E  4-  Cn  G]  cos(p7rs/n)  (50) 

(Mr,/a)„.,  =  D  l€iH  A  4-  «i9  C  4-  «2o  E  4-  €21  G]  cos(p7rs/n) 

(V»a),|.*  =  D  [€22  A  4*  *23  C  4"  *24  E  4"  *2.''i  G]  sin(p7rs/n) 

Again,  no  details  are  given. 

It  is  apparent  that  there  are  four  boundary  conditions  to  satisfy,  rather  than  three  as 
before.  These  are: 

for  the  clamped  edge 

(w/a),,..  =  0  ;  (lAr)q..  =  0  ;  (^.)q.,  =  0  ;  (a^'r.).,,.  =  0  (51) 

for  the  simply -supported  edge 

(w/a),,..  =  0  (^.)q,.  =  0  ;  (M,)q..  =  0  ;  (Mri(/fl)q,«  —  0  (52) 


Full  details  will  be  supplied  upon  request. 
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and  for  the  free  edge 

(Mr).,.,  =  0  :  (M,),,.,  =  0  ;  (Mr,/a)„.,  =  0  ;  (V-a),,.,  =  0,_  (53) 

For  each  of  these  cases,  a  4  X  4  frequency  determinant  or  transcendental  equation  is  obtained  and 
the  zero  values  of  tliis  determinant,  found  by  iteration,  correspond  to  the  desired  eigenvalues. 


4.0  RESULTS 

Twenty  eigenvalues  were  found  for  each  of  the  three  problems,  five  for  each  set  of  modes 
with  one,  two,  three,  and  four  half  waves  normal  to  the  simple  supports.  Of  each  five,  three  were 
symmetric  anti  ■  other  two  were  anti-symmetric  in  the  direction  parallel  to  the  simply-supported 
edges. 


For  the  non-conforming,  12-degree-of-freedom  element,  n,  the  number  of  elements  on  a 
side  was  varied  from  2  to  20  in  steps  of  2.  Some  larger  values  of  n  were  chosen  in  special  cases.  These 
results  are  given  in  Tables  5,  6,  and  7.  In  a  few  places,  especially  for  small  values  of  n,  it  was  found 
impossible  to  find  eigenvalues  because  the  characteristic  equations  had  imaginary  roots.  Blanks 
in  the  Tables  indicate  these  points.  It  was  also  impossible  to  find  solutions  for  the  higher  modes  for 
small  values  of  n.  No  great  effort  was  made  to  find  these  points  by  other  means,  since  the  results 
would  be  very  inaccurate  in  any  case. 

Since  exact  solutions  fcr  these  problems  are  available,  plots  of  error  versus  the  number  of 
elements  are  given  in  Figures  5,  6,  and  7.  Negative  error,  that  is,  solutions  below  the  exact  values, 
are  plotted  with  a  solid  line,  while  positive  errors  are  given  by  dashed  lines.  It  may  be  seen  that 
most  of  the  error  plots  tend  to  asymptote  to  a  straight  line  as  n  becomes  large.  The  slope  of  the 
straight  line  is  the  rate  of  convergence,  and  for  this  element  is  —2.  In  a  few  cases,  it  can  be  seen 
that  this  slope  has  not  been  reached,  even  though  the  number  of  elements  on  a  side  is  very  large. 

For  the  clamped  and  simply-supported  cases,  the  non-dimensional  parameters  sometimes 
start  above  the  exact  solutions,  cross,  and  then  converge  from  below.  This  cross-over  point  occurs 
at  small  values  of  n.  Similar  behaviour  is  observed  in  the  free  case,  but  the  cros-s-over  point  some¬ 
times  occurs  at  relatively  large  values  of  n.  Furthermore,  for  eight  modes  such  cross-overs  have  not 
occurred,  even  though  n’s  as  large  as 40  are  used.  It  is  not  known  whether  or  not  these  cross-overs 
will  ever  occur  for  these  modes. 

Percentage  errors  in  the  solutions  found  for  10  elements  on  a  side  are  given  in  Table  8. 
With  one  exception,  the  errors  are  all  greater  than  1%,  and  often  over  10%,  even  though  a  large 
number  of  elements  are  used. 

Little  difference  was  found  between  the  three-root  and  the  two-root  solutions  for  the 
clamped  and  simply-supported  cases.  Some  values  changed  in  the  fifth  or  sixth  place.  For  the  free- 
edge  case,  the  results  were  significantly  different  and  the  two  solutions  are  compared  in  Table  9 
for  some  of  the  eigenvalues.  As  expected,  the  eigenvalues  are  nearly  the  same  at  large  values  of  n. 
The  three-root  solutions  are  somewhat  better  at  low  values  of  n  when  the  solution  is  converging 
from  below,  but  are  worse  when  the  solutions  are  converging  from  above. 
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As  mentioned  previously,  Walz,  Fulton  and  Cyrus'"’  have  shown  that  the  rate  of  con¬ 
vergence  for  this  element  should  bo  n  “  for  the  simply-supported  case.  They  have  also  given  an 
expression  foi  the  error  in  the  solution  for  large  values  of  n.  Figure  8  gives  a  comparison  of  the 
Walz,  et  al.  predictions  of  error  with  the  actual  errors  found  for  a  few  typical  modes.  It  can  be 
seen  that  their  error  estimate  is  quite  good  for  large  values  of  n  (n  =  20)  for  the  lower  modes,  but 
not  as  good  for  the  higher  ones.  Nor  are  their  predictions  accurate  for  smaller  values  of  n. 


It  may  be  concluded  that  the  convergence  rate,  n-"^,  predicted  by  Walz,  Fulton  and  Cyrus'"’ 
for  the  simply-supported  case,  is  verified  herein,  and  that  this  convergence  rate  apparently  applies 
to  the  other  two  boundary  conditions  as  well  (except  for  some  modes  in  the  free  two  sides  case).  It 
must  be  emphasized,  however,  that  this  value  only  holds  for  very  large  values  of  n  and  does  not 
apply  for  array  sizes  that  are  readily  useable  in  the  direct  stiffness  method.  For  these  realistic  arrays, 
it  is  clearly  impossible  to  say  anything  general  about  convergence  rates. 


The  results  for  the  three  cases  found  using  the  16-degree-of-freedom  conforming  plate 
element  are  given  in  Tables  10,  11,  and  12.  Again,  20  eigenvalues  are  presented  for  each  case  and  the 
number  of  elements  on  a  aide,  n,  was  varied  from  2  to  10.  Solutions  for  small  values  of  n  were  more 
difficult  to  obtain  for  this  element.  Firstly,  the  characteristic  equation  often  had  imaginary  roots 
for  small  values  of  n,  and  this  prevented  solutions  from  being  obtained.  Secondly,  the  method 
broke  down  whenever  the  number  of  half  waves  in  the  simply-supported  direction  equalled  the 
number  of  elements  on  a  side.  To  complete  these  Tables,  eigenvalues  for  the  2  X  2,  3  X  3  and  4x4 
arrays  were  found  using  the  direct  stiffness  method.  These  results  are  indicated  by  a  star  in  the 
Tables.  In  many  places  results  were  found  using  both  methods,  and  the  values  always  agreed 
perfectly. 


Plots  of  absolute  error  vereus  the  number  of  elements  on  a  side  are  given  in  Figures  9,  10, 
and  11.  It  is  interesting  to  see  that  all  these  curves  rapidly  approach  an  asymptote  as  n  increases. 
The  slope  of  these  asymptotes  is  —4,  which  is  the  same  as  the  rate  of  convergence  predicted  for  the 
potential  energy  in  static  problems. 

It  is  significant  that  there  are  slope  discontinuities  in  these  error  plots  for  small  values  of  n. 
This  indicates  that  great  care  must  be  taken  in  extrapolating  results  calculated  for  small  values  of  n. 


It  is  particularly  surprising  to  see  that  the  convergence  is  not  monotonic  for  several  modes  in  the 
free  two  sides  case.  Since  this  is  a  conforming  element,  monotonicity  is  expected.  Closer  examination 
of  the  requirements  for  monotonic  convergence  indicates  that  finer  arrays  of  element  must  be 
contained  within  coarser  arrays  to  ensure  monotonicity;  i.e.,  results  from  2X2  and  4x4  arrays 
should  converge  monotonically,  but  results  from  2X2  and  3x3  arrays  would  not  necessarily 
have  to  converge  monotonically. 

Percentage  errors  in  the  solutions  found  for  10  elements  on  a  side  are  given  in  Table  8. 
These  errors  are  less  than  1%  for  all  but  one  mode  and  are  especially  small  for  the  lower  modes.  In 
general,  they  are  all  one  order  of  magnitude  or  morc  smaller  (up  to  three  orders  of  magnitude  in 
some  cases)  than  the  errors  found  using  the  non-conforming  elements.  This  Table  clearly  demon¬ 
strates  the  superiority  of  the  conforming  element. 
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5.0  CONCLUSIONS 

A  closed  form  type  of  analysis  using  shifting  E-operators  has  been  developed  to  study  the 
dynamic  convergence  of  finite  plate-bending  elements.  Two  rectangular  elements  have  been  studied, 
the  12-degree-of-freedom  non-conforming  element  and  the  16-degree-of-freedom  conforming 
element.  Three  dynamic  problems  involving  all  three  types  of  boundary  conditions  have  been  studied. 
Exact  solutions  were  found  for  these  problems,  so  that  convergence  of  the  finite  element  solutions 
could  be  considered.  Twenty  eigenvalues  were  considered  for  each  of  the  problems. 

The  closed  form  analysis  worked  well,  and  large  arrays  of  elements  could  be  studied  with 
little  computational  effort.  It  was  found  that,  with  few  exceptions,  the  non-conforming  element 
solutions  converged  from  below  the  exact  answers  for  large  values  of  n,  the  number  of  elements  on  a 
side,  at  a  rate  of  n  However,  the  array  size  needed  for  such  convergence  was  at  least  20  X  20  or 
larger,  so  this  convergence  rate  does  not  apply  for  array  sizes  that  are  readily  useable  in  the  direct 
stiffness  method.  The  convergence  rate  of  n”"  predicted  by  Walz,  Fulton  and  Cyrus  ’  for  square 
plates  simply  supported  all  round  was  confirmed,  but  their  estimates  of  the  error  magnitudes  were 
not  good  for  small  values  of  n. 

The  conforming  element  solutions  were  found  to  converge  to  the  exact  solutions  from 
above  at  a  rate  proportional  to  n~^  for  values  of  n  larger  than  6.  This  is  the  same  rate  as  predicted 
for  the  convergence  of  potential  energy  in  static  problems.  Slope  discontinuities  in  the  error  plots 
were  found  for  small  values  of  n,  indicating  that  great  care  must  be  exercised  in  attempting  to 
improve  calculated  results  by  extrapolation, 

A  comparison  of  the  errors  involved  in  using  these  two  elements  showed  that  the  conform¬ 
ing  element  was  far  superior  to  the  non-conforming  clement  in  both  magnitude  of  error  and  rate  of 
convergence. 
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FIG.7:  ERRORS  IN  NON- CONFORMING  SOLUTIONS  FOR  PLATE 
FREE  TWO  SIDES 
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FIG.9:  ERRORS  IN  CONFORMING  SOLUTIONS  FOR  PLATE  CLAMPED  TWO  SIDES 
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APPENDIX  A 

FINITE  DIFFERENCE  SHIFTING  E-OPERATORS 

(*••  R«f*r«nc«  IS) 

In  one-dimensional  space,  r,  these  operators  are  defined  as 

EJ  w,  =-  w,+i, 

and  provide  a  concise  notation  for  the  function,  w,,  specified  at  finite  points.  In  two-dimensional 
space,  r  and  s,  the  operators  are 

Ef  E„  w, ,  “  + 

These  operators  obey  the  rule 

F,  (E,'  'E.;  s'-'  1>™  =  a‘'  b™  F,  F.  ib"! 

I'sing  this  rule,  finite  difference  or  finite  element  equations  may  be  traasformed  into  ordinary  expo¬ 
nential  equations.  For  example, 

iEr‘  -E^siniTTr  H)  -  (E;'  -  EM-ite‘""  -  e"*' /2  ) 

_i  j  .f-"  '<  _  e”  e'"  "  -  le"  ^  -  e“*  '*)  e'”'  )  2 

-  -  2  sin  ITT  Hi  cr«;)rr  K> 

A  list  of  the  transformations  needed  in  the  analysis  are  given  below. 

SHIFTING  E-OPERATOR  TRANSFORMATIONS  k  it  any  integer 


(e;" 

-  E‘.  e“‘'  -  - 

2  sinh  k<T  e"'* 

iEp-‘ 

-  E‘  e"*’  = 

2  cosh  kcr  e"*’ 

(E;'‘ 

-  Efi!  e”  =  - 

2  i  sin  kiT  e"' 

-f  E^  i  e'"  = 

2  cos  k<r  e"'’ 

—  Ep'  sinh  ffp 

=  —  2  sinh  k(j  cosh  ffp 

-r-  Ep!  sinh  <Tp 

-  2  cosh  k<7  sinh  op 

fEp-' 

-  Ep  .'  cosh  up 

-  —2  sinh  ko  sinh  op 

(Ep-'' 

f  Ep  )  cosh  ap 

=  2  cosh  ko  cosh  op 

'.Ep^ 

-  Ep)  sin  <rp 

=  —  2  sin  ko  cos  op 

(E,:' 

+  Ep  .  sin  ap 

=  2  cos  ko  sin  op 
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(Ep '■  -  Ep)  cos  ffp  =  2 

(Ep'‘  +  Ep)  cos  ffp  =  2 

(Ep'‘  —  Ep)  (  —  1)”  sinh  <tp 
(Ep'‘  +  Ep)  (  -  I)'*  sinh  ap 
(Ep'‘  -  E^)  (-l)‘*cosh  cp 
(Ep‘‘  +  Ep)  (-l)’’cosh  ap 
E;-''  e*'”  = 

Ep  cos  ap  -  cos(p-(-  k)<r 


sin  k<r  sin  trp 
cos  k<T  cos  ap 

-  2  sinh  ka  (  —  1)’’  cosh  ap 

=  —  2  cosh  kff  ( -  I)**  sinh  ap 
=  2  sinh  ka  (  —  I)**  sinh  ap 

=  —  2  cosh  ka  (  - 1)’’  cosh  ap 
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APPENDIX  B 

DEFINITIONS  FOR  X  AND  e 

The  constant  terms  appearing  in  the  expressions  for  (^r)r.«  and  (^B)r,i.  (t*q.  (36)  and  (37)) 

are  as  follows: 

Xi  =  —  sinh  «  (  [  39  cos(p7r/n)  4-96-1-7  (116  cos(p7r/n)  4  274)  ]  [cosh  «  (17  cos(p7r/n) 

4  22)  4  28  cos(pir/n)  4  68  4  10  7(cosh  «  4  2)  (3  cos(p7r/n)  —  4)] 

4  28  7  sin-(p7r/n)  [39  cosh  «  4  96  4  7  (116  cosh  a  4  274)]}  /X? 

Xo  =  sin  0  {[39  cos(p7r/n)  4  96  4  7  1116  cos(pir/n)  4  274)  ]  [cos  fi  (17  cos(p7r/n)  4  22) 

4  28  cos(p7r/n)  4  68  4  10  7  (cos  3  4  2)  13  cos(p7r/n)  -  4)]  4  28  7  sin-(p7r/n) 

[39  cos  4  96  4  7(116  cos  8  4  274)  ]  }  /X« 

X;i  is  the  same  as  Xi  with  sinh  «  replaced  hv  —  siiih  J  and  cosh  «  by  —  cosh  P,.  The  denominator 

is  now  Xn. 

\4  -  —  sinlp7r/n)  [|39  cosh  «  4  96  4  7  (116  cosh  n  4  274)]  [cos(p7r/n)  (17  cosh  n  4  22) 

4  28  cosh  a  4  68  4  10  7  (coslpr/n)  4  2)  (3  cosh  a  —  4)  ] 

-  28  7  sinh^o  (39  coa(p7r/n)  4  96  4  -ydlO  cos(px/n)  4  274)  ]}  /X7 

Xr,  =  —  sin(p7r/n)  (  [39  cos  4  96  4  7  1116  cos  3  4  274)  ]  [coslpTr/n)  117  cos  3  4  22) 

4  28  cos  3  4  68  4  10  7  lcoslp7r/n)  4  2)  13  cos  3  —  4)  ] 

4  28  7  sin^tJ  [39  coslp7r/n)  4  96  4  7III6  coslp7r/n)  4  274)  ]  )  /Xs 

Xr.  is  the  same  as  X.i  with  sinh  «  replaced  by  —  sinh  p  and  cosh  n  by  —  cosh  $.  The  denominator 

is  now  X;). 

X7  =  -  784  7- sinh'-'a  sin'-^lpTr/n)  —[cosh  a  117  coslp7r/n)  4  28)  4  22  coslp7r/n)  4  68 
4  10  7  IcoslpTT/'n)  4  2)  13  cosh  a  —  4)}  (cosh  n  117  coslp7r/n)  4  22) 

4  28  coslpTT/n)  4  68  4  10  7  Icosh  «  4  2)  13  coslp7r/n)  —  4)) 

Xj,  =  784  7-  sin'-'/J  sin*lp7r/n)  ~  (cos  8  117  coslp7r/n)  4  28)  4  22  cos(p7r/n)  4  68 
4  10  7  IcoslpTT/n)  4  2)  13  cos  3  -  4)]  (cos  3  117  coslp7r/n)  4  22) 

4  28  Cf)slp7r/n)  4  68  4  10  -ylcos  3  4  2)  13  coslp7r/n)  -  4)} 

Xn  is  the  same  as  X7  with  sinh  n  replaced  by  —  sinh  p  and  cosh  n  by  —  cosh  p. 
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The  constant  terms  appearing  in  the  expressions  for  (Mr),,,H,  (Mh),,,b  and  (V*a),,.s  (eq,  (40) 
are  as  follows: 

<1  =  X|  |Zi  sinh  rt(q  -1)4-  7j>  sinh  «q  +  10  >  Z.-j  (3  sinh  «(q  —1)  —  4  sinh  aq)  ) 

+  14  7  (2  cosh  «(q  -1)  +  3  cosh  wq)  sin(prr/n) 

4-  Z4  cosh  «fq  -1)  —  Zr,  cosh  «q  4-  7  [Za  cosh  «(q  — 1)  4-  Z7  cosh  aq  ] 

(■'  =  X2  |Z,  sin  0(q-l)  +  Zo  sin  (3q  4-  10  7  Z,i  (3  sin  3(q-l)  -  4  sin  0q)  ] 

4-  14  7  Xr,  (2  cos  /3(q  -1)  4-3  cos  ;3q)  sin(pn-/n) 

-f  Z4  cos  0(q  -1)  -  Zr,  cos  (3q  4-  7  (Zr,  cos  /3(q  -1)  -f-  Z7  cos  /3ql 

«n  =  Xa  [Z,  (-1)"-'  sinh  ?(q-l)  -f  Zo  (-1)*'  sinh  £q  -f  10  7  Z3  (3  (-l)s-'  sinh  {(q-1) 

-  4  (-!)'•  sinh  iq)]  4-  14  7  Xr.  (2  (-!)•'-'  cosh  f(q-l)  4-  3  (-1)''  cosh  fq)  sin(p7r/n) 

4  Z4  (-1)"-'  cosh  ^(q-1)  -  Zr.  (-1)"  cosh  ^q  4-  7  [Zr  (-1)''-'  cosh  {(q-1) 

4-  Z7  (-1)'’  cosh  fql 

<4  =  -  14  Xi  7  [sinfp7r/n)  (2  sinh  «(q  -1)  -  3  sinh  aq)  ] 

+  X4  [Zr  cosh  oCq  — 1)  -f  Z<i  cosh  aq  -f  10  7  Zio  (cosh  o(q  — 1)  -f  2  cosh  aq)  ) 

-  sin(p7r/n)  [39  cosh  a(q  -1)  4-  96  cosh  aq  -  7  (116  cosh  a(q  -l^)  -f  274  cosh  aq)  ] 

«r,  =  -  14  X2  7  [sintpTT/n)  (2  sin  d(q  -1)  —  3  sin  0q)  ] 

4-  Xr,  [Zr  cos  0(q  -1)4-  Za  cos  ^q  4  10  7  Zio  (cos  0(q  -1)4-2  cos  0q)  ] 

-  sin(p7r/n)  [39  cos  0(q  —1)  4-96  cos  0q  —  y  (116  cos  8(q  —1)  -I-  274  cos  ^q)  ] 

fr,  =  -  14X3  7  [sin(p7r/n)  (2  (-1)''-'  sinh  $(q-l)  -  3'(-l)''  sinh  tq)  ] 

4-  Xfi  [Zr  (-!)''-•  cosh  ?(q-l)  4-  Zo  (-1)''  cosh  £q  -f  IO7  Z,o  (  (-1)"-'  cosh  $(q-l) 

4-  2  ( -  I)**  cosh  Sq)  1  -  sin(p7r/n)  [39  (-!)'•-'  cosh  $(q  -1)  -f  96  ( -1)^  cosh  ^q 

-  7  (116  (-1)"-'  cosh  5(q-l)  -f  274  (-!)•>  cosh  ^q)  ] 

«7  =  -  X(  (Z4  sinh  a(q  -1)  -f  Zr.  sinh  aq  4  7  (Z«  sinh  a(q  -1)  -  Z7  sinh  aq)  [ 

-  X4  8in(p7r/n)  [39  cosh  a(q  -1)  4-  96  cosh  aq  -f  7  (116  cosh  a(q  -1)  4-  274  cosh  aq)  [ 

-  3  Z2  cosh  atq  —1)  —  Zti  cosh  aq  —  7  (Z12  cosh  o(q  —1)4  Z13  cosh  aq) 
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<8  =  -  Xa  IZ4  sin  /i(q  -1)  +  Z5  sin  /Sq  +  7  (Zr,  sin  ^}(q  -1)  -  Z7  sin  /3q)  1 

—  Xr,  sin(p7r/n)  [39  cos  /3(q  —  1)  4-  96  cos  /Jq  +  7  (116  cos  |3(q  —  1)  +  274  cos  /3q)  1 

-  3  Za  cos  /3(q-l)  -  Zn  cos  (3q  —  7  (Zia  cos  /?(q  -1)  +  Zn  cos  0q) 

<9  =  -  Xn  IZ4  (-!)■>-’  sinh  $(q-l)  4-  Zr,  (-1)’'  sinh  ^q  4-  7  (Zr,  (-1)""^  sinh  $(q-l) 

-  Z7  (-!)’’  sinh  fq)  ]  -  Xr,  sin(p7r/n)  [39  (-1)''-'  cosh  $(q-l) 

4-  96  (-!)'»  cosh  £q  4-  7  (116  (-!)■'-'  cosh  £(q-l)  4-  274  (-1)"  cosh  ?q)  ] 

-  3  Za  (-1)''-*  cosh  ^(q-1)  -  Z,i  (-1)'’  cosh  {q  -  7  IZia  (-1)'’“'  cosh  {(q-1) 

4-  Zi3  (-!)'»  cosh  tq] 

where 

Zi  =  17  cos(p7r/n)  4-  28  ;  Za  =  22  cos(pir/n)  4-  68  ;  Z;,  =  cos(p7r/n)  4-  2  ; 

Z4  =  39  cos(p7r/n)  4-  96  ;  Zr,  =  24  cos(p7r/n)  4-  111  ;  Zg  =  116  cos(p7r/n)  4-  274 

Z7  =  199  cos(p7r/n)  4-  461  ;  Zh  =  17  cos(p7r/n)  4-  22  ;  Z9  =  28  cos(p7r/n)  4-  68  ; 

Zio  =  3  cos(p7r/n)  -  4  ;  Zu  =  204  cos(p7r/n)  -  474  ; 

Zia  =  394  cos(p7r/n)  4-  1226  ;  Z13  =  1226  cos(p7r/n)  4-  3454 


